#pragma once

#include <vector>
#include <map>
#include <cmath>
#include <algorithm>

// permutation entropy

namespace entropys {

template<typename T> std::vector<int> argsort(const std::vector<T> &array)
{
    const int array_len(array.size());
    std::vector<int> array_index(array_len, 0);
    for (int i = 0; i < array_len; ++i)
        array_index[i] = i;
    std::sort(array_index.begin(), array_index.end(),
        [&array](int pos1, int pos2)
        {
            return (array[pos1] >= 0 && array[pos2] >= 0) ? (array[pos1] < array[pos2]) : (array[pos1] > array[pos2]);
        });
    return array_index;
}

inline int factorial(int n)
{
    long long res = 1;
    for (int i = 2; i <= n; i++) {
        res *= i;
    }
    return res;
}

inline double pe(std::vector<double> x, int m, int t)
{
    using namespace std;
    if (t > m)
        t = m;
    vector<vector<double> > X;
    int length;
    if (t == 1)
        length = x.size() - m + 1;
    else
        length = int((x.size() - m + 1) / t) + 1;
    for (int i = 0; i < length; i++) {
        X.push_back(vector<double>());
        for (int j = 0; j < m; j++)
            if (i * t + j <= x.size() - 1)
                X[i].push_back(x[i * t + j]);
            else
                X[i].push_back(-1);
    }
    vector<vector<int> > index;
    for (int i = 0; i < X.size(); i++) {
        index.push_back(vector<int>());
        vector<int> ar = argsort(X[i]);
        for (int j = 0; j < X[0].size(); j++)
            index[i].push_back(ar[j]);
    }
    map<vector<int>, int> symbol_map;
    for (int i = 0; i < index.size(); i++) {
        if (symbol_map.count(index[i]) == 0) {
            symbol_map.insert(pair<vector<int>, int>(index[i], 1));
        }
        else {
            symbol_map[index[i]] ++;
        }
    }
    double result = 0;
    for (auto iter = symbol_map.begin(); iter != symbol_map.end(); ++iter) {
        double p = 1.0 * iter->second / (1.0 * (x.size() - (m - 1) * t));
        result -= p * log(p);
    }
    result = result / log(1.0 * factorial(m));
    return result;
}

inline double mpe(std::vector<double> x, int tau, int m, int t)
{
    int y_len = x.size() - tau + 1;
        std::vector<double> y(y_len);
    for (int j = 0; j < y_len; j++) {
        double sum = 0;
        for (int i = j; i < j + tau; i++) {
            sum += x[i];
        }
        y[j] = sum / tau;
    }
    return pe(y, m, t);
}

}
